d c = 4 is the upper critical dimension for the Bak-Sneppen model 
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Numerical results are presented indicating d c = 4 as the upper critical dimension for the Bak- 
Sneppen evolution model. This finding agrees with previous theoretical arguments, but contra- 
dicts a recent Letter [Phys. Rev. Lett. 80, 5746-5749 (1998)] that placed d c as high as d — 8. 
In particular, we find that avalanches are compact for all dimensions d < 4, and are fractal 
for d > 4. Under those conditions, scaling arguments predict a d c = 4, where hyperscaling 
relations hold for d < 4. Other properties of avalanches, studied for 1 < d < 6, corroborate 
this result. To this end, an improved numerical algorithm is presented that is based on the 
equivalent branching process. PACS number(s): 64.60.Ak, 05.40.-a, 05.65.+b. 



The Bak-Sneppen evolution model jjj has received con- 
siderable attention as an archetype of self-organized crit- 
icality Q , which has been put forward as a general mech- 
anism leading to many non-equilibrium scaling phenom- 
ena observed in Nature S. Although the Bak-Sneppen 
model was originally developed as an attempt to inter- 
pret paleontological data indicating co-evolutionary ac- 
tivity in biological evolution jj], Q] , it has also been used 
to interpret power law distributions in quiescent periods 
between earth-quakes Q. Its generic mechanism of ex- 
tremal dynamics || has even inspired an algorithm to 
approximate combinatorial optimization problems Q. 

In a recent Letter || , numerical results for several crit- 
ical exponents of the Bak-Sneppen model in high dimen- 
sions were presented. On the basis of that study, and cer- 
tain scaling arguments, Ref. || concluded that the upper 
critical dimension, where those exponents obtain mean 
field values, was d — 8. Avalanches, cascading through 
the system via nearest-neighbor activation events, form 
domains that were found to be fractal for d > 2. These 
two findings are surprising: On one hand, avalanches 
would proceed on a filamentary domain where surface 
sites outnumber bulk sites. On the other hand, activ- 
ity would have to return mostly — by pure chance — to 
sites already within that domain to keep it from growing 
linearly with time, as required in the mean-field limit. 

The claims for Ref. Q are in stark contrast to earlier 
derived scaling relations (see Tab. 1 in Ref. ||), which 
would predict an upper critical dimension of d — 4. For 
instance, the scaling relation for the avalanche cut-off is 
a = l-r + d/D. Its mean-field value is a = 1/2 [§ |o|, 
and the avalanche distribution exponent attains r = 
3/2 §,0, while Ref. explicitly derived that D = 4 
is the mean field value for the avalanche dimension expo- 
nent, implicating d — 4 as the upper critical dimension. 
But the scaling relations in were derived under the 
clearly stated assumption that avalanches would always 
be compact for d < d c , and activity within them would 
proceed homogeneously over their domain, with many 



returns to each site. Ref. || argues that avalanches are 
ramified already for d > 2, allowing for a d c > 4. 

In this letter we test the claims of Ref. || by indepen- 
dent means. To this end, we have developed an alterna- 
tive algorithm, which has the added benefit of improv- 
ing temporal cut-offs by more than two decades while 
eliminating finite lattice size effects entirely. As a result, 
we find the assumptions underlying the scaling theory 
in |1 intact. In particular, our numerical simulation re- 
sults show that avalanches in d < 4 are compact, and 
that scaling exponents, when asymptotic behavior is ex- 
tracted, take on mean-field values for d = 5 and 6, very 
much in contradiction to Ref. Q . Thus we show that the 
Bak-Sneppen model possesses an upper critical dimen- 
sion, d c = 4, below which hyperscaling relations hold, 
similar to equilibrium critical phenomena. In fact, we 
have studied a whole range of different properties of the 
Bak-Sneppen model for 1 < d < 6 pi] ]. Here, we focus 
directly on those properties relevant with regard to the 
upper critical dimension. 

The Bak-Sneppen model is very easy to state It 
consists of random numbers A r between and 1, occu- 
pying sites r on a (i-dimensional lattice. At each update 
step, s, the smallest random number \ m in(s) is located. 
That site as well as its 2d nearest neighbors each receive 
new random numbers drawn independently from a flat 
distribution on the unit interval. This update step is re- 
peated at the site with the next X m in{s + 1), and so on. 
The process inevitably evolves toward a self-organized 
critical state in which almost all A r are larger than a 
critical threshold A c (see Tab. |), while those A r below 
are part of an avalanche of activity. The current mini- 
mum X m in{s) is the "active" site while those A < A c are 
"unstable" and potentially active; all A > A c are "sta- 
ble". These avalanches are critical and their properties 
have been described in terms of scaling relations 0. 

The obvious way to implement the model is to specify 
a lattice, place a random number A r on each site r, keep 
an efficient list of those numbers, ordered by size, and re- 
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peatedly replace the smallest \ m in(s) and its neighbors. 
In this method, knowledge of A c is a priori not necessary 
to determine the critical properties of the avalanches. 
This is, in fact, the algorithm employed in Ref. M and 
some other, low-dimensional studies. But in higher di- 
mensions feasible length scales in the pre-defined lattice, 
with a fixed number of sites N — L d , exponentially dimin- 
ish with dimension. For instance, Ref. [|| used N = 2 16 , 
giving L — 1Q in d — A (only L = 4 in d=8), or a maximal 
spatial separation of r max = dL/2 = 32, and a temporal 
cut-off of maximally s co ~ r^ ax « 10 6 , since D < 4. 

A more efficient way to describe an avalanche utilizes a 
generalized branching process |12, ijjj that has been shown 
to be exactly equivalent to the Bak-Sneppen model. In 
the update procedure, the values of the new random num- 
bers that enter the system are unrelated to each other, 
their predecessors, or any other number in the system. 
Thus, any number that has no chance of becoming ac- 
tive in itself (i. e. those with A > A c ), will not affect the 
dynamics of the critical avalanche through its particular 
value. Now, assume that we are at the beginning of a 
critical avalanche, i. e. all numbers in the system are 
above A c . All we need to know to describe the ensuing 
critical avalanche is the existence of a smallest number 
in the system to initiate the avalanche. In updating that 
number and its neighbors, either no number below A c 
will be created and that critical avalanche happens to 
be empty, or it will produce numbers below A c and will 
progress through more updates until no more numbers 
below A c remain at some later time step. Thereafter, 
a new critical avalanche will start at a new, unrelated 
location. To describe an individual avalanche, we can 
assume that it started at the origin, and at all times we 
need only keep track of those numbers below A c . Lat- 
tice addresses appear only as parameters characterizing 
exclusively those numbers currently unstable. While the 
number of sites covered by an avalanche n cov can grow 
as much as linear in time s, currently unstable sites at 
most grow like s 1 / 2 [||. Thus, our temporal cut-off due 
to memory constraints N behaves like s co ~ iV 2 . 

Since our results depend so strongly on results from 
numerical simulations, we describe our improved method 
in great detail. To simulate this process we create a list 
C of currently unstable numbers, C — {(A r ,r)|A r < A c }. 
We initialize C with a single entry (A r = 0,r = 0). At 
each update we first remove the smallest A r in C and any 
of its neighbors in the lattice, if those happen to be in 
C. Then, we draw a new number for each of those sites, 
but only sort into C those numbers below A c by storing 
(A r , r). Addresses r could extend arbitrarily far from the 
origin, thus eliminating any spatial cut-off. It is very easy 
to keep C compact and sorted according to A. In fact, it 
is sufficient to sort in a new A r linearly from the bottom 
(instead of using a heap, say): the dynamics quickly leads 
to a list in which almost all numbers are very densely 
packed just below A c , and almost every number inserted 
at the bottom moves only a few steps through C. 

Without explicit lattice reference, it is not easy for this 



TABLE I: List of the critical threshold values A c , and 
the exponents D and /i, for 1 < d < 6. The relative error 
for A c is of the order of 10 -5 . The errors in the exponent 
are about 1%. The values for d = 1,2 are taken from 
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D 



/' 



1 0.66702 2.43 0.41 

2 0.328855 2.92 0.68 

3 0.201190 3.35 0.905 

4 0.142412 4 (log-corr.) 0.99 

5 0.110020 4 1 

6 0.089752 4 1 



algorithm to check the minimum's neighbors in the lattice 
which, if unstable, would need to be removed from C. To 
be sure, we would have to search C at each update to 
eliminate those neighbors, which would be unreasonably 
time consuming. Instead, we utilize a procedure similar 
to hash tables |k| : When a lattice address r with A r < A c 
is stored in some entry Ck, r is mapped into an index 
i = i(r) for a large, sparse array A, \A\ ~> \£\. Ai in 
turn stores the index k of Ck, or is empty otherwise. 
When activity returns to the neighborhood of r, z(r) is 
calculated and Ai can be checked at once to track an 
unstable number in C. It is crucial that i(r) is unique, of 
course, but since |.4| is finite there exist lattice addresses 
r ^ r' such that i(r) — i(r'). Those conflicts are rare 
between any two currently unstable sites, if \A\ 3> \C\. 
They can be resolved by moving to the next available 
index i(r) + I, where / is the offset to the nearest free 
entry in A (for important details, see Ref. pj]] ). 

We have first used this simulated branching process 
to extrapolate for the values of A c ]ll| given in Tab. | 
by simulating avalanches at various A below A c (see 
Ref. ||). This extrapolation uses the domain covered 
by an avalanche, n cov , requiring C to record all currently 
and previously unstable sites. At A c , n CO v can increase as 
much as linear with s, quickly exhausting memory (see 
below). Fortunately, for A well below A c , avalanche du- 
ration and coverage are quickly cut off. 

We have run up to 2 x 10 6 critical avalanches in each 
d, 1 < d < 6, to determine their statistical properties. 
We used |£| max = 2 16 and \A\ = 2 20 , easily sufficient 
to run avalanches up to s max = 10 8 -C |£| 2 lax update 
steps. [Avalanches were terminated at s max due to time 
constraints and due to the error in A c , see Tab. |^.] A is 
already much larger than the lattice used in Ref. [pj , but 
we point out that even on a fixed lattice with N = 2 24 
sites in d = 6, the largest distance is r max = 48, i. e. s co < 
10 6 . In contrast, our algorithm has s co ~ l^lmax ~ 5 x 10 9 
in any dimension, yielding significant spatial correlations 
for r > 100 even in d = 6 (see Fig. |l|). 

In a long critical avalanche, almost all unstable num- 
bers cluster densely in L with values just below A c , 
which puts high demands on a random number gener- 
ator. Typically, in an ongoing avalanche that has grown 
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FIG. 1: Plot of s/(r) D as a function of (r) for the 
distribution of avalanche activity in 1 < d < 6. For d < 4 
we used the fitted values Dd=i = 2.43, Dd=2 = 2.92, and 
Dd=3 = 3.35. For d — 5 and 6 we used the mean- field 
value D — 4. For d — 4 we used D — 4 as well, but also 
a factor of ln(r) as indicated by Eq. ([!]). 

to 10 4 unstable numbers at some time s, 10 3 of those 
are placed within AA/A C « 10~ 4 just below A c . Using 
a = 1/2, those numbers can expect to survive up to 
s ~ (AA/A C ) _1 / CT ~ 10 8 update steps. To record tempo- 
ral correlations accurately, it is crucial to resolve these 
numbers accurately. This requires numbers that are suf- 
ficiently random for more than 8 digits! Our data has 
been obtained with a sophisticated 64-bit random num- 
ber generator provided to us by the authors of Ref. Q . 

First, we discuss the data for the distribution of 
avalanche activity, P(r,s), which leads to the avalanche 
dimension exponent D. To find P(r, s), we record the 
instances of having activity at a (positive) distance r rel- 
ative to the origin at update step s. Similar to a random 
walk, the moments of the distribution define the expo- 
nent D via (r q ) s ~ s q / D . In Fig. [I] we have plotted the 
data as s/ (r) D vs. (r) for d ^ 4. For d < 4 we have used 
the value of our best fit for the exponent D || |ll| as 
given for each dimension in Tab. |. For d > 4, the fitted 
value corresponded to the mean-field result ju|, D = 4. 
For d = 4, we conjecture a simple asymptotic form for 
the scaling behavior, 

(r) 4 

(* = 4). (D 

Thus, in Fig. [IJ unlike for the other dimensions d, for 
d = 4 we have plotted sln(r)/(r} 4 . For all dimensions 
the data levels out horizontally for increasing (r). Loga- 
rithmic factors, as seen in Eq. ([!]) for d = 4, are common 
for scaling behavior at the upper critical dimension [ jT5| . 

Similar evidence for d c — 4 is provided by the data 
for the backward avalanche exponent t^ 11 , which is re- 
lated to the avalanche distribution exponent through 
t = 3 — r^ n ||. This particular scaling relation does 
not rely on avalanches being compact and is valid both 
above and below the upper critical dimension. The data 
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FIG. 2: Plot of P£ u (s)/s~ 3 / 2 as a function of l/ln(s) for 
backward avalanches for d — 1 to 6 and for the mean-field 
multi-trait model. For d < 4 the asymptotic behavior, 
l/ln(s) — > 0, is characterized by a rapid approach to 
zero, while for d > 4 the ratio does not approach zero. In 
the marginal case, d = 4, P^ u (s)/s -3 / 2 appears to vanish 
as some power of 1/ ln(s). (The multi-trait data has been 
shifted upward by 0.1 for clarity.) 

exhibits mean-field behavior, t^ 11 = 3/2(= r), already for 
d = 5 and 6. In Fig. || we show the data for the backward- 
avalanche distribution reduced by the mean field scaling 
behavior, P£ ll (s)/s~ 3 / 2 as a function of 1/ ln(s), for d = 1 
through d = 6, and for the multi-trait model fltj|| , which 
analytic calculations show exhibits mean field behavior. 
The curves for d < 4 clearly approach zero rapidly for in- 
creasing times s on this scale, indicating that r^ n > 3/2 
in those dimensions. The corresponding data for d > 4 is 
clearly bounded away from zero, comparable to the data 
from the multi-trait model, where r^ n = 3/2 exactly. 
This data also shows that there are significant correc- 
tions to scaling and cut-off effects, even for the multi- 
trait model. In this case, simply fitting a power law over 
the largest region available from the simulation data can 
give misleading results due to systematic curvature on 
increasing scales. (Fitted values for r^ 11 are discussed in 
Ref. [0.) In d = 4, the data is again consistent with 
some logarithmic correction to scaling behavior. 

These findings receive strong theoretical support. The 
only assumption in the scaling theory of Ref. § is that 
the avalanches are compact below the upper critical di- 
mension where each site that is visited in an avalanche is 
typically visited many times. Our numerical results in- 
dicate that avalanches indeed remain compact for d < 4. 
To test this explicitly, we have determined the probabil- 
ity distribution P(n cov , R) of having a finished avalanche 
that has covered a domain n cov with a radius of gyra- 
tion R. As mentioned above, n cov may grow linear with 
s, rapidly exhausting memory for C. Thus, we used 
smax = |£Ux = 2 20 « 10 6 and \A\ = 2 22 . The ex- 
ponent d cov is defined through the moments of the dis- 
tribution P(n COV7 R) via (n cov ) ~ R d c°v f or large R, pro- 
viding a measure of the fractal structure of avalanches if 
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FIG. 3: Plot of (n cov )/R d as a function of l/\n(R d ) 
obtained from the coverage n cov and the gyration radius 
R of each avalanche in d = 1 to 6. Here, the quantity 
appears to converge to a finite limit as l/ln(i? d ) — » 
for d < 4, signaling the compactness of avalanches. In 
turn, for d > 4 it rapid approaches to zero. In the 
marginal case, d — 4, (n cov )/R d appears to vanish lin- 
early in 1/ \n(R d ). 
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FIG. 4: Plot of the coverage, reduced by its mean-field 
scaling, (n cov )/s as a function of s. The data for d = 5 
and 6 is horizontally level over many decades, consistent 
with mean-field behavior. The data for d = 3 exhibits 
scaling consistent with its value of D. The deviations 
from mean-field scaling in d = 4 are insignificant. 
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d cov < d. In turn, avalanches are compact if d cov = d. 
Using our numerical results, we plot in Fig. [| the quan- 
tity (n cov ) / R d as a function of l/ln(i? d ). If avalanches 
are fractal (dcov < d), we would find that asymptotically 
this quantity should approach zero (like #-(<*-*=<>«)). i n 
fact, we find that this quantity clearly remains finite for 
large R for all dimensions d < 4. For d > 4, we find a 
rapid approach to zero on this logarithmic scale, indicat- 
ing fractal behavior. In d = 4 this quantity appears to 
vanish as well, but as a linear function of 1/ \n(R d ). Sim- 
ilar to Eq. (jl]), we find that the numerical simulations 
indicate marginal behavior in d = 4 with logarithmic 
corrections. Thus, in d = 4 avalanches are marginally 
compact (a "fat fractal"), and the conditions underlying 
the scaling theory in Ref. || are upheld for all d < 4, 
implying d c — 4. 

We demonstrate the consistency of our argument by 
also measuring the scaling of the coverage with the life- 
time of an avalanche, (n col ,)~s M ||. Since there should 
be only one characteristic length for a compact avalanche, 
we expect that (r) ~ R, and thus, n = d/D for d < D, 
and /x = 1 in the mean-field limit. In Fig. we plot 
the coverage reduced by its mean-field scaling, (n cov )/s, 
as a function of duration s. Clearly, the data for d = 5 
and 6 is in perfect agreement with mean-field behavior, 
while for d = A the deviation from mean-field behavior are 
minute, ^ moas ps0.99, even without considering the effect 
of logarithmic factors. For d = 3 we find /i~ 0.905, quite 
consistent with the value predicted from fi = d/Dm 0.896 
using our measured value D = 3.35 in d — 3 (see Tab. |). 
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